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I. INTRODUCTION 


In recent papers Karmarkar [Refs. 1,2] has presented a new 
method for the solution of linear programming (LP) problems: His new 
solution technique moves from some feasible starting point across the 
interior region of a polytope that is defined by the problem constraints. 
He shows that the number of steps to find an optimal solution with his 
technique is polynomially bounded. 

In contrast, the simplex algorithm which, is widely used for the 
solution of LP problems finds an optimal solution by moving from vertex 
to vertex on the polytope. This is known, in the worst case, to require 
an exponential number of steps [Ref. 3]. 

The polynomial bound of the projective algorithm makes the new 
solution technique very appealing to researchers. The theory of the 
new algorithm seems to be widely accepted among experts, while 
Karmarkar's claim that his algorithm is 50-100 times faster than the 
simplex method has met skepticism [Ref. 4]. 

In this study a variant of the projective method is implemented, 


and some well known test problems are solved. 
A. THE BASIC PROJECTIVE METHOD 


Folowing Karmarkar [Ref. 1: p. 4], the number of steps of the 
algorithm depends on R/r, where R is the radius of the sphere 
circumscribing the polytope, and r the radius of the inscribed sphere. 


Assume a general LP of the form 


Min ely 
Sy tea 


UI 
oO 


(LP1) 


* 
IV 
oO 


—_ 


With the assumption that the sum of its variables has an upper bound, 


and with the proper scaling of variables, a convexity constraint 


itx =1 (1.1) 
can be added to LP1. This transformation maps the LP onto a unit 
simplex S$ whose center is at eee cra, L/ Te Then, a transforma- 
tion is performed such that a feasible interior starting point is mapped 
onto a). 

If Bidet) is the largest sphere with center a, that can be 


inscribed into the simplex S, and “B(a,,R) is the smallest sphere 


circumscribing S, then 
Rye =n. (122) 


By restricting a solution x to remain inside the largest inscribed sphere 
Exerc), the method achieves in one iteration a reduction in the differ- 
ence between the current objective value and the optimal objective value 
by a factor of (1 - 1/n). Following Shanno [Ref. 5] a simple proof is: 
Let 


f" = min clx, Seeeoe AX = oD, (eure 

f = min ol x, Kee iar), (20'S =! 16 itx = Ts ea) 

f =mineclx, xe BiGe Kee x = b, lx = 1. (alae 
Then, 

cede tS ean a eta. = f, (iG) 


and with equation 1.2 


Sh a T 


oO Ci eee Ce a.) ele) 
From that 
(£- £ )/(eta, - f°) < (1 - I/n). (1.8) 


If a linear objective function could be maintained at each iteration, 
it follows that an upper bound on the number of steps required to find 
an optimal solution is of O(n). Unfortunately, the projective transfor- 
mations needed to continue the algorithm result in a nonlinear objective 


funetion: 
B. DERIVATION OF THE ALGORITHM 
1. The Canonical Form 


Suppose we have a linear programming problem of the form 
LP1. Karmarkar's method requires that this LP be transformed into the 


following canonical form, 


Min el x 

s.t. Ax = Q (LP2) 
ilx = 
x Ve U0 


where the optimal solution value is zero. With the assumption that the 
sum of the variables of an LP is bounded above and subsequent scaling 
by this bound, a convexity constraint (equation 1.1) can be added to 
LP1. In practice this can be a problem because choosing too large an 
upper bound may cause numerical problems. 

LP2 requires that the nonhomogeneous system of equations 
Ax = b be made homogeneous. Karmarkar [Ref. 1: p.34] proposes a 

T 


transformation that would transform the i-th equation 2 Sa b; to 


ay = bx; = Oo (i) 
This transformation has the disadvantage that when b is dense the 
sparsity of A will be lost. Karmarkar requires the optimal solution 
value to LP2 to be zero together with a special stopping criterion 
(equation 1.26), to prove the polynomial bound. To achieve the zero 


objective value the optimal objective function value f of the 


10 


untransformed LP has to be known in advance, and the following trans- 


formation made 
oly - f= ely - Tx = (ce -f1)Tx . (1.10) 


Karmarkar concludes [Ref. 2: D. 387] that if the minimal objective value 
determined by the projective algorithm is not equal to zero, the original 
problem must be either infeasible or unbounded. Another method for 
making the transformation from LPl to LP2 is discussed in 
Chapter f1.5B. 


Ze The Projective Transformation 


Let x9 > 0 be a known feasible solution to LP2. The invertible 


projective transformation 


Dlx Gein 
y= —— 
18 p-lx 
maps any xX, such that ifx = 1 and x2 0, onto y such .that 
ily = 1 and jy 2 0. D is a diagonal matrix whose entries are 
ese... xX, )- 


The transformation maps the LP2 unit simplex in x-space onto 


another unit simplex in y-space. The point x° is mapped onto 


eee t/n 10 the center of the unit simplex in y-space. The inverse of 
the transformation (equation 1.11) is given by 
Dy CipalZ> 


itp VY 


After the transformation, LP2 can now be restated as 
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Caw oy 
Min 
itp y 
s.t. Av Diya — ae (LP3) 
ily Ss 
y 2 0 


Strictly speaking LP3 is not a linear program because its objective 
function is a rational function of y. Karmarkar [Ref. 1: page 18ff] 
shows that it is sufficient to consider a linear approximation to the 


objective function of LP3, which leads to 


Min ec!D y 
ee NID (LP4) 
ily = 
y 2 0. 


3. Optimization Over a Sphere 


A solution to LP4 is now restricted to le within a sphere with 


center at y® and radius ar, where 
r= 1/(n(n-1)) 1/2 . (1.13) 


r is the radius of the largest sphere that can be inscribed into the unit 
simplex, and qa is a constant such that O<q<l. q provides a margin 
which ensures that the algorithm doesn't select a point outside the 
sphere due to round-off error. By convexity, an optimal solution will 
occur at the boundary of the sphere. Thus, the additional constraint 
can be stated as an equality. Now we have_the following version of the 
LPs 


1 


Min clp y 


s.t. ADYy = QO (LPS) 
ily = 1 
(y-y®) F¢y-y®) = a ?r2 


Before the problem can be solved, one more transformation 
that moves the center of the sphere to the origin is useful. Let 


y= y+ y®, and eliminate the constant terms Ay — 0. iTy° = 1 and 
Typ,,0 
D 








ec Dy~. For convenience, define B by 
Aa 1.14 
om [AD (1.14) 
a 
and the gradient ¢ of the objective function of LP5 by 
Esc D . (1.15) 
Then LP5 can be restated as, 
Min @ly 
Se) a = QO (LP6) 
yl¥y S oe 


In order to solve LP6, note that an initial feasible solution to 
LP6 is Y= 0 (which corresponds to y = y® in LP5). This is also the 
TS a a“r@ 


obtained by finding a direction of maximum rate of ascent @ that is 


center of the sphere y An optimal solution to LP6 can be 
feasible with respect to By = 0, and moving in direction -€ (maximum 
rate of descent) a distance ar from y= 0 to the boundary of the 
sphere. 

A feasible direction of maximum ascent is found by orthogo- 
nally projecting ¢ onto the null space of B (see [Ref. 1: page 17] ), 


i.e., the following problem has to be solved in terms of Cy 
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Min Go = gaUeuermemmec 


pp 
=p 


Sie Be, = i. 


Define the Lagrangian L(c,,, A) Soy 


~ alsa _ oxl 
L(c,,,A) = & 2c Cp 


+ 


p 


Cc 


p 


. om 2 
Min (Ie - chil 2) 


p 


First-order optimality conditions for the Lagrangian yield 


20 cams BIAx = 0, 


and 


After multiplying by B and dropping 2Bec 


-2B@ = BBA 


p 


=O (equation 1.19) 


Assuming (BBL)-1 exists, we can solve for A 


A = -2(BB1) 1 Bz. 


Substituting equation 1.21 into equation 1.18 gives 


Cs ech Bt (BB!) 1 Be. 


The direction of maximum rate of ascent is then 


Cc = Cp / Il Cyl : 


Thus, the optimal solutions to LP6 and LPS are 


Y = - Are, 
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(1.16) 


(1.17) 


(1.18) 


(CIES ibs) 


we get 


(1320) 


C2 1p 


Gir 22) 


(ieoy 


2) 


Ve = yy -are Cie 25) 


respectively. 


Finding Cy is the key part of Karmarkar's method because it 
involves the major portion of the computational work. Solving equation 
ie22 or Cy can be viewed as solving a linear least-squares problem 
(see Chapter II.C.). 

With the optimal solution to LPS in y-space, the method 
proceeds by transforming that solution back into x-space by use of 
equation 1.12. Next it defines a new matrix D _ and iterates until it 
reaches the stopping criterion [Ref. 1: p. 14] 

(el xkys(elx°) < 279 Gee) 
where q is a termination parameter. Note that ol x = 0 at optimality. 


Table 1 shows an outline of Karmarkar's basic method. 
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TABLE 1 
ALGORITHM 1 


Input Problem: SiZe@ ae eee n 
Coetfiiesz-ent. Ma tai ee A 
Cost FUnNC@Ul Oni ee c 
Initial Feasible Soluticnm: 3 
Termination Parameter....... q 
Feasibility Parameter....... ‘e,4 
Begin k = 1 
5, ae 
f° = oly? 
f = ow 
r = 1/(n(n-1))71/4 
While =) (f/f =e a 
= diag (x) 
(D7txy/(itp-tx) = 1 


Ol << O 
II 


U9) 
II 





fe) 
I 


¢ - Bl ppl)-tpe 
2 Cp/ I Cyl 

=A = lies 

(Dy) /(1T Dy) 

= ely 

= kt#il 


A Hh KM <M OD 
I 


End(While) 
End 
Output SO LUC i.OM. 5 slau a eee eee x 
Objective Function Value...f£ 


l[teration Counts... eee oe k 
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IT. A VARIANT OF THE PROJECTIVE METHOD 





A. INTRODUCTION 


The following sections outline a variant of the projective method 
that was proposed by Wood [Ref. 6]. It has a practical method of 
bringing an LP into the canonical form, and uses the non-linear objec- 
tive function of LP3 to find a gradient c' which corresponds to ec of 
equation 1.16. Similar to the simplex method, the proposed algorithm 
uses a ratio test to determine a feasible step length. 

Other practical features of the proposed method include the relax- 
ation of the requirement to know the optimal objective value in advance, 
and exploitation of sparsity of A. The implementation is especially 
concerned with controlling fill-in during the solution of the normal 


equations. 
B. DERIVATION OF THE MODIFIED ALGORITHM 


Given an LP problem of the form LP1l, we use a single artificial 


variable x,,4,, to attain initial feasibility: 


Min (CoG ep 
s.t. Bece (AED a = 8 (LP7) 
SF ne | 2a Oe. 


where M is the cost for the artificial variable and x° is the initial solu- 


tion with 
eer sh 
car aap = le Onl) 


The transformation has added one more column to the problem. 
To get a homogeneous right-hand side in LP? we introduce an 
additional variable, and apply the following projective transformation, 


whose inverse is given by 
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(XX a 7) = 2) Ca a oe (22) 

C2 ae (n+2)(X,%y4q)/(17 (Xp a) 4D) Oe 

Ki = (n2)) Geka ase (2.4) 
Let the artificial column be denoted by a, i.e. a= -(Ax°-b). The 


transformed problem can now be restated, with the exception of the 


objective function, in Karmarkar's canonical form, 


aE t t t 
M,0O)(X',xX x ) 
, ee y) p) ) nie | n+? 
Min ae 
nt+2 
St. (A,a,-b)(xX',x .41,% ,39 (LP8) 
Hi (pace x aan) = nt2 
Ce a) = 0. 


The above transformation adds yet another column to the problem, but 
has the advantage that it doesn't change the sparsity of A, as opposed 
to Karmarkar's proposal. Rather than projecting the LP problem onto a 
unit simplex, it is projected onto an (n+2)-simplex. This is done to 
improve numerical stability. 


The projective transformation, equation 1.12, is applied to LP8 


which gives 


Min) (C+ 4M 0) (G sya ee) one a 


CULM ID vivian) a=) © (LP9) 
LAGE ye eee = M2 
(Ye) nse) 2 0 < 


The gradient (compare with equation 1.16) of the objective function of 
LP9Y, evaluated at the initial feasible starting point fA = 11 is propor: 
tional to 
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Be (cjdy)-2- c,d Md. .,,-c x). (2.5) 


The step of normalizing Cp in Algorithm 1 is replaced by a ratio 
test. Let 


jmax = argmax{(C),);}. (225) 
This allows the algorithm to make a step outside the inscribed sphere, 
but maintains feasibility by restricting a solution to lie inside the 


simplex. Then, the update of y becomes 

v= 1 - PCy/ (Cp) imax : (22) 
where p is a parameter to maintain feasibility such that O<p<1. Table 2 
Shows the modified algorithm. 


C. THE LINEAR LEAST-SQUARES PROBLEM 


Computation of the projected gradient c, during every iteration 


Pp 
accounts for most of the computational workload in any algorithm based 
on Karmarkar's method. Solving equation 1.22 can be viewed as solving 


the following linear least-squares problem, | 
Min (Ile! - BAM)? (2.8) 
where Bt is an (n+2)xm matrix with ms<(nt2) assumed. 
If rank(B1)=m, then the solution to (2.8) is given by the solution 
to the system of normal equations 


BBIA = Be’. (2.9) 


The projected gradient Cp is then the residual vector of the least- 


squares problem (2.8), i.e., e 


c - c! = Br ; (2.10) 


Ihe, 


TABLE 2 
ALGORITHM 2 


Input Problem S22 © et. eee n 
Coeffietent) Mare 1c ere A 
Cost Function gee eee ot 
Initial Feasible Solution...x® 
Termination Parameter....... 6) 
Feasibility Paranever sia 

Begin 
£° = Oo 


While (£95) £ (eee) 





D = diag (xX) 
fe == f(x) 
y = (n+2)(D7+x)/(itp" tx) 
' = 
Cc Vyt(y) 
A D 
B = 
ab 
es = ¢'- BlpBl)-tse' 
1 eal argmax(Cy)j 
y = Ls Pc,/ (Cy) imax 
x See ee Oye ye 
k = eee 
End (While) 
End 
Output 50.1 UWE Onl: 2 eee 2 ee x 
Objective Function Value...f 
LteraCion: Counme.. ..ee ee k 
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BB! is symmetric and positive definite, given that Bl has full column 
rank. 

As Heath [Ref. 7: p. 499] points out, the ideal choice for solving 
the normal equations, given full rank, is Cholesky factorization. If Bt 
is not of full column rank the cross-product matrix will be singular and 
(BBI)71 ceases to exist. BB! could become nearly singular if BI is 
near rank degenerate. Shanno [Ref. 5: p. 25] shows that this will 
happen if the optimal solution is degenerate, i.e. as the optimum is 
approached numerical problems arise and Cholesky factorization is likely 
to fail. 

One problem that cannot be avoided when solving the normal equa- 
tions is the fact that the P-condition number of BB! is the square of 
that of Bt [Ref. 8: p. 223], so that when BI is already ill-conditioned 
it may be impossible to find an accurate solution to equation 2.9. 

Another important consideration when computing BB! is that the 
sparsity in Bl will not automatically guarantee sparsity in BB!. In 
fact, the addition of variables in LP? and LP8 has added two possibly 
dense bottom rows into Br, seus, BBI will be completely dense. 
However, one can cope with that by initially omitting the dense rows in 
Bt from the computation of BBL, and later updating the solution to 
equation 2.9 using procedures similar to the ones described in 


ieee o: Pp.) 9s-65]. 
D. IMPLEMENTATION 


Algorithm 2 has been implemented in FORTRAN H_ (Extended) 
Opt(2) on an IBM 3033 AP under VM/CMS. Al floating point arithmetic 
is performed in double precision. The program is designed to accept 
different solution modules from available software packages for solving 
the linear least-squares problem. 

Input data sets are in standard MPS format. The numerical values 
of the non-zero elements of the constraint matrix A are stored column- 


wise in a real array. For versatility a full set of pointers are defined: 


Le column index 


Ze: R row index 
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3. AP location of first non-zero element in a column 
RP location of last non-zero element ina row 
5. LINK location of next non-zero element (backwards) 


in a wow 


Brameller, Allan and Hamam [Ref. 10: p. 104-110] give’ a comprehensive 
discussion of sparse storage schemes. 

The transformation of an LP problem into the canonical form adds 
one dense row and two dense columns to the A matrix. The dense row 
originates from the convexity constraint equation 1.1, the first dense 
column stems from the single artificial variable that is needed to attain 
initial feasibility, and the second dense column is added to make the 
system of equations homogeneous (see LP8). The artificial column is 
updated with the residual of the current solution as long as the total 
infeasibility is above a specified threshold. 

As mentioned earlier, dense rows in Br yield BBL completely 
dense. With the following method, this problem can be alleviated. The 
method applies to any number of dense rows in Bt, but in this study 
we are only concerned about two dense rows, namely the ones that 
result from the transformations that are performed to get from LP2, to 
LP8 via LPT. Consider the projection problem of equation 1.16 in the 


folowing form, 


Mingle semen lee Gail} 
Sb Be, = 0 

Replace Cp by Zz, and let 
B = | 5 (2. 
c= (cance (2513) 
i = (21, 250 : (2.14) 
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where By is an (mxnj), Bo an (MxN9) matrix, and No is the number of 


dense columns in B. Equation 2.11 now becomes 
Cie Cie Zoll 5)" * (Ne', - 2,11 5)7 (2.15) 
2 2" 2 il 1" 2 
Sint By24 = ~BoZo 


and when separated, 2.15 becomes 


2 : i Z 
ff Cle’ - zoll 9)? + Min (Ihe'y - 2,Il9) (2.16) 
Min 21 


Ls - 
St: By,2, = ~BoZo . 


Solving the inner minimization first, the Kuhn-Tucker conditions yield 


Cy 2 Z1 = BAA ) (Zane) 
and 
ByZ; = =B5Z5 (2.18) 


Multiplying 2.17 with B, we get, 
= T 

Bycy cs By2Z1 me B,B, A ° (2.19) 
Substituting 2.18 into 2.19 gives 

B,c', + Boz) = BB IA. (2.20) 
Let Ag solve 

B,c', = B,B,!A (221) 

Pe ie 

and let Aj solve - 


(Bg), = 3 (8). a ea ea Ce 
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where (Bo); is the j-th column of Bo, then 2.22 corresponds to solving 


the matrix equation system 
By BEES (2523) 
2 1 age 


where A is an (mxny) matrix. Then, the solution to A for the inner 


minimization is, 
9 + Ads . as 
Substituting 2.24 into 2.17 gives 

ce!) = 2, = Ba ore are | (2.25) 


To simplify notation, let h = Bie Ag and H = B TA 


Then equation 2.25 becomes 
Z1 = cy ~ h = HZ5 e (ZezZ0)} 
Substituting 2.26 into 2.15, the overall minimization then becomes 
As ete Z 2 
Min (le'p ~ Zyl 9)? + (Ih + Hzgil 9) (2.27) 
Ze 


which is a simple least-squares problem. The first-order Kuhn-Tucker 


optimality conditions yield, 

c's - hiH = (I + HtH)z, (2.28) 
which gives, solving for Zo, 

Zo = (1 Helee(c', ah JHee (2.29) 
With the solutions to 2.29 and 2.26 we have the desired result, 


Cy = 2% = (Ze ; (2230) 
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In practice the following procedure is followed: 


le: compute B,B,7 

Sa Laccor Beer eeoemusing Cholesky factorization 

ee solve B,B,TA = uO ie les fee Ubibniest Og Mabie) WAV 

4. compute h = BTA, 

Bis solve BB 1A See yee —solution is ia, 

.. solve B,B, 1A Sao yee ter SOlUCTCmE 1s a) 

Tes compute H = BTA 

ex compute (Cle Hly)-t, (note that H'H is a 2x2 
matrix if there are 2 dense rows in Bl) 

9. Gomipace (C5) hiH) 


10. compute Z>5 
ii Ccomoute Z4 


The given procedure is efficient since the factorization of BB! is 
computed only once and the same system is solved three times using 
different right hand-sides each time. 

As a stopping criterion for the algorithm the following rule is 


used, 


IF argmax (i : ce D = | hor (Zc) 
j 


where t is a real constant. The convergence criterion 


CHI Cy M/\otx|) cant (232) 


mentioned by Lustig [Ref. 4: p. 12] can also be used to terminate the 


algorithm. 


III. SOLVING THE LINEAR LEAST-SQUARES PROBLEM 
A. INTRODUCTION 


In this study, the linear least-squares problem is solved by 
explicitly computing and solving the normal equations. Although the 
normal equation approach experiences problems when applied to il- 
conditioned or near rank deficient matrices, it behaves acceptably with 
sparse and well-conditioned matrices [Ref. 11]. 

The numerical methods for solving normal equations generally fall 
into two classes, direct and iterative methods. One representative of 
each class is considered in this study. Little can be said as to which 
class of methods is better, except that direct methods are more attrac- 
tive in terms of computational work, and iterative methods may require 


less storage [Ref. 12: p. 11]. 
B. CHOLESKY FACTORIZATION 


The method implemented is given in George and Liu [Ref. 12], and 
uses Cholesky factorization with a minimum-degree ordering to solve a 
large sparse positive definite system of equations. Since BBL is only 
guaranteed to be positive semidefinite, a modification to the Cholesky 
factorization algorithm is considered in Chapter IV.A.3. to accommodate 
the semidefinite case. 

The minimum-degree algorithm is a reordering heuristic which 
attempts to reduce fill-in during the factorization phase. The reordering 
phase is entirely symbolic; it amounts to a symmetric row and column 
permutation of BB! which corresponds to reordering the columns in 
Bt. During this phase, BB! doesn't have to be computed numerically; 
only its structure has to determined. Also, the factorization is first 
performed symbolically, thus allowing a static data structure for the 
Cholesky factor L. An outline of the phases of the algorithm is given 
in Table 3. See also Heath [Ref. 7: p. 499]. 
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. TABLE 3 
MINIMUM-DEGREE ORDERING / CHOLESKY FACTORIZATION 
ALGORITHM (MDOC) 
Determine the nonzero structure of BB!. 
Paci ma permutation Matrix =P such that PBB pt 
has a sparse lower triangular Cholesky factor 
Diy: 


Factor PBBip? symbolically and set up the data 


Secrelure Lor i, 
Compute pBBipt 
Factor PBBIp! = uy! numerically. 


pemve  Uze— PBD §) (back substitution). 
Ab 
L 


numerically. 


y=2Zz (forward substitution) . 





1. The Minimum-Degree Ordering Heuristic 


The heuristic finds an ordering of a symmetric matrix such 
that fill-in is low when the matrix is being factored. The basic idea is, 
at each (simulated) factorization step, to permute the part of BBL 
remaining to be factored so that a column with the fewest nonzeros is in 
the pivot position. The implementation consists of six subroutines that 
are given in George and Liu [Ref. 12: pp. 124-137]. The subroutines 
accept as input the adjacency graph associated with BBL represented 
by an adjacency structure, and return as output a symmetric permuta- 
tion of BBL given as a permutation vector for the columns of Br. 

Let G=(X,E) be the adjacency graph of BBL, where X is 
the set of nodes, and E is the set of edges. Then, the nodes corre- 
spond to the variables of the least-squares problem, i.e. the columns of 
Bl. Two nodes x and y are said to be adjacent if {x,y} is an edge in 


E. The adjacent set of Y, Y¢X is defined and denoted by 


ATi 


Adi(Y) = (% € X-YI(X, YieE forisomes, eae (oe) 


An eaten list for xXeX is a listporalif@nedes” meee {<))- 
Finally, an adjacency structure for the graph G is the set of adjacency 
lists for all xeX [Ref. 12: pp. 37-41]. The particular adjacency 
structure used in the ordering heuristic stores elements in each adja- 
cency list in contiguous locations. An entry point array to the first 
element in each list allows access to the list. 

The ordering heuristic is based on graph theory, and involves 
the notion of elimination graphs, quotient graphs, reachable sets and 
indistinguishable nodes. George and Liu [Ref. 12: pp. 92-124] may be 


consulted for more details. 
oe Factorization and Solution 


The components of the lower triangular Cholesky factor L of 
BB! are computed using the so called "inner product form" algorithm 


[Ref. 12: p. 20]. The elements of L are given by, 


j-1 

= ia = 2 1/2 j= 

Ly = (ey 7 Pp for j=1,2,...,m (3.2) 
ee 

lij = (ei; - 2 rene EE fOr t=) 17) +2 een (353) 


where the Ci; are the elements of BBI. 
After the factorization has been computed, the following two 


linear systems have to be solved (see also Table 3): 


Lz = Pb! , (3.4) 
and 
Lily =z. : | (3.5) 


28 


Solving system 3.4 by back substitution involves the use of "inner 
prediwets” [Ref. 12: p. 25], defined by 


ee noe om 
Zi = (b; 2, Liz) /h; for t= 1s 2 ' (3.6) 


where bY stands for the i-th element of the right-hand side of (3.4). 
For this case L must be accessed row-by-row. 

System 3.5 is solved by forward substitution using "outer 
products" [Ref. 12: p. 26], defined by 


wae Z,/1,; hOomot=i2,...,m (3.7) 
as: 2) ae ee 2 a Yiiey p> -otn Dd . 
For the latter case, LI must be accessed LOWY COW mor  scoluml-= py — 


column instead. 
C. INCOMPLETE CHOLESKY FACTORIZATION 


This method is an implementation by Ajiz and Jennings [Ref. 13] of 
the incomplete Cholesky conjugate gradient algorithm (ICCG), whose 
theory is given in [Refs. 14,15]. The algorithm requires that the 
coefficient matrix of a set of simultaneous linear equations be symmetric 
and positive definite. It consists of two distinct parts, one being the 
Cholesky factorization, which can be complete or incomplete, and the 
other a conjugate gradient iteration to solve a preconditioned linear 
system, where the Cholesky factor serves as the preconditioner. 

Golub and Van Loan [Ref. 16: pp. 373-377] point out that precon- 
ditioning is essential for obtaining good convergence rates with conju- 
gate gradient methods. The convergence rate is closely linked to the 
P-condition number, which is the ratio of the largest to the smallest 
eigenvalue. It was mentioned earlier that the condition number of BBt 
will be the square of the one of Br. IH-conditioned problems have 
large condition numbers, and hence slow convergence. Preconditioning 
is a process of transforming a linear system so that its P-condition 


number is improved [Ref. 17: p. 979]. 


TRS, 


Consider again the original linear least-squares problem (2.9) with 


a generic right-hand side b' and x the unknowns, 
BB! x = b'. (3.8) 


Then, equation 3.8 can be preconditioned with a transformation matrix 


L giving 

Lipptin tls = Lt, (3.9) 
On 

Lipp in ly = 5b, (3.10) 
where y = Lix and b = L 1b’. According to Ajiz and Jennings 


[Ref. 13: p. 950] the ideal choice of transformation matrix L is the 
Cholesky factor of BBI, since 


Lippitt =1, @kit) 


provided one could perform exact arithmetic. 

The objectives of the incomplete factorization phase of the algo- 
rithm are to transform BB! as close as possible to I, and to reduce 
fill-in in the factor L. This is accomplished by discarding some off- 
diagonal coefficients during the factorization, whose magnitudes fall 
below a preset threshold limit. The result of this operation is an incom- 


plete Cholesky factor L that must satisfy 

BB! =aph er (Beai2) 
where C is the matrix of elements omitted from the factorization. 
Unfortunately, omission of elements from the factorization. process can 


destroy the positive definiteness property and hence lead to a break- 


down of the process. Ajiz and Jennings [Ref. 13: pp. 950-951] have 
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shown that introducing diagonal modifications in C will retain the posi- 
tive definiteness property. The matrix C will be symmetric and will 
have diagonal elements that are greater than or equal to zero. Thus, C 
is a positive semidefinite matrix. Assuming BBL to be positive definite, 
adding C will result in LLE also being positive definite. 

Convergence of conjugate gradient iterations is only guaranteed 
for the positive definite case. Thus, a modification to adapt the 
Cholesky factorization to the positive semidefinite case as with the 
direct method may cause slow convergence. The modification is consid- 
ered in Section IV.A.3. Table 4 gives an outline of the ICCG 


algorithm. 


TABLE 4 
INCOMPLETE CHOLESKY - CONJUGATE GRADIENT (ICCG) 
ALGORITHM ; 
Obtain L, an incomplete Cholesky factor of BBL. 
Solve Lbw=eb'  fommbs byetorwardwcubstitution. 
Solve L~tppiy-ly = b for y by conjugate 


gradient iteration. 


Determine xX by back substitution in Lix = va 





1. The Incomplete Factorization 


The procedure presented here is given in Ajiz and Jennings 
[Ref. 13: pp. 951-952], and Jennings and Malik [Ref. 14: pp. 
310-313]. To see how the elements in column j of L are computed 
consider the following. From matrix equation 3.12 a typical elemental 
equation may be written as 
j-1 


perc cn > Sil (ied) (Geis) 


Ij 1) 1 Kel 


J) 
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where the subscripts for elements | refer to their positions in the lower 


triangular factor L, and the ej; are the elements of BB!. Let 


ste j-1 


“.. =e. - mals. 3 314 
Sei ij ey Lindi. . ( ) 


Then, assuming the elements of L have been computed for columns 1 to 
(j-1), all the ej 
where all the ei pass the rejection test, i.e., are to be retained. 


in column j can be computed. First consider the case 


Hence, by setting ¢4j-0 in equation 3.13, we get 


Lg = ea /]y - (3.15) 


In case any ei is to be rejected, lj, is set to zero, implying Gre iy: 


The off-diagonal term Cij implies diagonal additions Ci and c;, to the 


matrix C in order to have the positive semidefinitness property. The 
matrix C doesn't have to be stored in memory; only the diagonal addi- 


tions Ci are of further interest. 


Before any li; can be computed, 1;; has to be determined. With 


j) a 
all the diagonal: additions Cai that resulted from rejections of ei -during 


the computation of columns 1 to (j-1), an expression for 1; becomes 


Wea Gli Ge OS a USN (3.16) 
JJ Kejel JJ 


where e ;: is defined by 


JJ 
ais j-1 j-1 7 
“= ee: CK). . 
oo > a Cas me ee 5 (3217) 


and os) is the diagonal addition to Ci resulting from deletion of ek 


and erg the diagonal addition to c;, resulting from deletion of ea. 


The rejection operation tests the magnitude of an element en 
in relation to the current values of the corresponding diagonal elements 


é.. and ei respectively, whose values are given by 


Suse ., + See (3.18) 
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where the index k refers to the rows for which rejections in column j 


= 
ee 


have had a bearing on aT and 


J-i1 
Ci, = ii ne ey) (3.19) 


In equation 3.19, k refers to the columns for which rejections in row i 


have had a bearing on ey 


An element e'; is rejected if 


] 
ae Biieui > 20) 


where y is the preset rejection parameter. A choice of y=0 will retain 
all elements, thus leading to a complete Cholesky factorization. A choice 
of y=l1 will cause all off-diagonal elements to be rejected. Ajiz and 
Jennings [Ref. 13: p. 952] recommend that the rejection parameter be 
in the range 0.01<y<0.2 for effective incomplete factorizations. 

The diagonal modifications in C which result from the rejection 


of an element en are given by 


eo iat ie. ./a.. 1/4 

ail le kj! (e9j7ek? ’ (3.21) 
and by 

= le 1(6.,/e,,.)2/- . | (anes 


With the successive application of equations 3.14 in conjunction with the 
relection operation, 3.17, 3.16 and 3.15, all elements in column j of L 


are determined. 
2. The Conjugate Gradient Iteration 


The conjugate gradient method of MHestenes and Stiefel 
[Ref. 18] is applied to matrix equation 3.10. It uses the following 


vectors, the letter k indicates the k-th iteration, 


a) p(k) Conjiugace gradient vector 


b) r(%) residual vector 
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Cc ) y(¥) Solution) yeocecs 
ay u'*) product ofet) 5BE- uso sandio a 
The initial values for k=0 are y(9) = 0 and p(9) = r(9) = b. The algo- 


rithm for one iteration is as follows 
ufK) = po lppli-Tp(*) 

y(K+1) = y(K) 4 appl) 

r(k+1) 2 pC) 2 ay (0) 


ce ( KL) y TROL) 7p CX) y TEC) 


By, 


p(k+l) = p(k+1) 4 Bi p() 


The first step in the above algorithm is obtained without computing the 


transformed matrix explicitly by the following three operations, 


Ae Liy(k) = p(k) (back substitution) 
oe w(K) = pply(®)  (pre-multiplication) 
So. bulk) = w(K) (forward substitution). 


The algorithm is terminated when 
ir(K) 7 WHI < tolerance (3.23) 


where b equals the starting residual (9) since y (9) was chosen to be 


Zero. 


34 


IV. MODIFICATIONS AND COMPUTATIONAL RESULTS 
A. ALGORITHMIC MODIFICATIONS 
1. Iterative Improvements 


A provision to improve the solutions to equations 2.21 and 
2.22 has been implemented, because their residuals are often unduly 


high. The residuals are computed as 
r= B,B,1A = b! ; (4.1) 


where A stands for Xo, Ay and Avy respectively and b' is generic for 
By Ronee pnmand (Bo)o. 
If Ir;l> e for some i, where ¢ is set to, say, 10° the 


following systems are solved for \' 
B,B,+ A‘ = -r. | | (4.2) 
An improved solution is then obtained by adding A and X! 
B,B,1( A+ A’) = rtb'+(-r) = b’. CoD 


The improvement can be repeated if the residuals of equation 4.3 are 
still found to be too high. ‘With this modification the direct method 


(MDOC) has become a semi-iterative method. 
2. Removal of the Artificial Column 


The update of the artificial column with the residuals of the 


current solution 
Ges) bx' 49 = 3 | (4.4) 


fet 


Pear?) - r= Ax 1) - b - x/x'_,.*) (4.5) 
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on each iteration has been augmented, such that the artificial variable 
can be deleted from the problem when the infeasibility becomes small. 
After the deletion, only one dense column remains in B, and the solu- 
tion of the least-squares problem is simplified. The motivation for this 
modification is to avoid some of the numerical problems that arise from 


variables that approach zero. 
3. Positive Semidefinite BB! 


It can be shown that if, and only if a diagonal element ever 
goes to zero in a Cholesky factorization, BBL is positive semidefinite, 
not positive definite. Furthermore, all elements below the zero diagonal 
element must also be zero and a nonunique solution to the normal equa- 
tions can be obtained by setting the corresponding A; to 0. This solu- 
tion can be obtained by replacing the zero diagonal element with any 
positive value and continuing with the factorization. Restated, if a 
diagonal element li in the partially computed Cholesky factor is less 
than or equal to io © set 1 i=1, and 1,,=0 for i=j+1,j*+Z,...,mM in eee 
tions 3.2 and 3.3 respectively. This procedure is also useful to deal 


with numerical problems which arise from degeneracy near optimality. 
4. Weighted Homogeneity Variable 


The initial solution (Yuya) = 1 used to begin the 
projective algorithm is arbitrary, and in some sense, the "homogeneity 
variable" ¥° 442 is fundamentally different than the other variables. 
Consequently, a modification has been made to allow weighting the 
starting solution Vow differently from the other ar I=1,.52, ne leew 


yo 4078. With 11 y%=n+2, the other y°; are set equal, 
¥, = (nt2=s)/ (ei) i eee (4.6) 


It is hoped that such a weighting might lead to lower iteration counts. 


~ 
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B. TEST PROBLEMS AND COMPUTATIONAL RESULTS 


At first, a small test problem (TEST1) was used to verify the 
correctness of the source code. The two algorithms were then tested on 
seven problems which have been also used for testing by Lustig 


(Ref. 4: p. 19]. Table 5 shows some of the characteristics of the test 


problems. 
TABLE 5 
TEST PROBLEMS 
Name Rows Columns Logicals Nonzeros Density 
eo. t 6 2 6 10 Oxses 
AFIRO oy 32 19 95 O-103 
ADLITTLE 56 oy 41 322 O. O86 
SHARE2B 96 79 SS 901 0.119 
ISRAEL 174 142 174 PASTAS. WO) 
BRANDY jis 249 54 2204 0.046 
E226 VARS: Boe 190 2518 0.041 
BANDM 305 472 0 Zoog OO is 
Table 6 summarizes the computational results. Alle Cr URS imes 


represent the time in seconds to set up the major part of -the data 
structure, solve the LP, and write out a few parameters on each itera- 
tion and the solution. Times to read in the data from MPS format are 
not included. 

The convergence criterion equation 2.31 is used with t=0.05 for all 
test problems but TEST1 and SHARE2B, where t=0.001 and t=0.01 
respectively. No solutions have been obtained for the problems AFIRO, 
BRANDY and BANDM: the algorithm will not converge to the optimum. 

The feasibility parameter p (equation 2.7) is set to 0.9995, giving 


the best overall performance of the algorithm. Tests with p=0.9999 and 


SIT 


p=0.99999 on the data sets TEST1 and SHARE2B, indicate that with 
these higher values the number of iterations necessary to attain feasi- 
bility is reduced, but the total number of iterations remains about the 
same. In addition, the accuracy of the solutions deteriorates. 

An intermediate solinon is declared feasible, i.e. the artificial 
variable is removed from the problem, as soon as X,4,M is less than 
10°}. The choice of 10°! is arbitrary, and depends on the value of the 
optimal solution. Using alternate values of 1.0 or 10°2 makes little 


difference in the number of iterations necessary to attain feasibility. 


TABLE 6 
COMPUTATIONAL RESULTS 


MDOC ICCG 


Problem Iterations Iterations 


feasible optimal feasible optimal 


TESTI 1 6 1 6 
ADLITTLE 13 Dik 13 mAh 
SHARE2B 6 wy 7 a7 
ISRAEL 19 98 19 89 
E226 8 4g 9 ag 





Factorization failures plagued both algorithms before the modifica- 
tion to accomodate a semidefinite BB! was implemented. These failures 
occured near the optimum, e.g. E226, or with the ICCG algorithm when 
setting the rejection parameter to a value greater than zero. After 
implementing the semidefinite modification these factorization problems 
have been cured, but the ICCG algorithm now shows very slow conver- 
gence. For example, a solution to SHARE2B is obtained only after 
112.36 CPU seconds with y=0.015 and all other parameters unchanged. 


Thus, because storage is not at a premium, only complete factorizations 
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are used. Satisfactory numerical results with the incomplete factoriza- 
tion are obtained only with the trivial test problem TEST1; computation 
times are always inferior. 

The minimum-degree ordering works very well in practice. 
Computational cost seems to be moderate, e.g., 0.25 seconds for 
SHARE2B, and 3.12 seconds for BANDM. Table 7 gives densities of 
BB! and the Cholesky factors with and without reordering. Substantial 
reductions in storage requirements were achieved when doing the 


reordering. 


TAB EEe? 
DENSITIES 


Re-ordered Not re-ordered 


Nonzeros Nonzeros Nonzeros 


Name in BB? Density in L_ Density in L Density 


TESel 20 
AFIRO 62 
ADLITTLE 384 
SHARE2B 271A 
ISRAEL 22 7 
BRANDY 2853 
E226 2823 
BANDM B72 


BZ 26 
ee A LOW? 
.241 411 
ES] TOZt 
ioe 11433 
mS 2 3429 
syle SHSISEe, 
. 080 4660 


Be vA) 
neo5 194 
23 816 
5 AES, ial 34 
o> | we 7s3 
ses 71760 
. 146 O/ 35 
100 SOS 


B52 
3 
Paz 
a2 o3 
eos 
oa! 
. 430 
.688 


@ © @ 0 0O GO © 
Gao OVO © [O80 ©€& 
GeO © @ OC O80 ©& 


The densities and number of nonzeros in L or BBL are relative 
to a symmetric half of a matrix; the diagonal elements are not 


included. 





Computational results indicate that iterative improvements are not 
always an absolute necessity. Residuals of a current solution tend to be 


high at the start of iterations, probably due to a less than optimal 
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choice of M. Doing a few (2 or 3) iterative improvements at this stage 
stabilizes the computations until the algorithm gets close to the optimal 
solution. Numerical problems near the optimum are so severe that even 
allowing a prohibitively high number like 20 iterative improvements has 
no apparent influence on the quality of the solution. 

The mildest form of numerical difficulty near the optimum is slow 
convergence, e.g. ISRAEL and E226. For SHAREZB, a high quality 
solution is obtained with no numerical difficulty. This data set is chosen 
to test the weighted homogeneity variable modification. S is given 
several values in the range 1 to 100. Iteration counts range from 20 to 
37. The quality of the solutions is sometimes reduced, however. Only 
S=39_ (27 iterations) and s=50 (28 iterations) give high quality solutions 
with low iteration counts. Values of s greater than 75 creat conver- 


gence failures. 
C. CONCLUSIONS 


The low iteration counts for some of the test problems are prom- 
ising, although CPU times seem to tell the difference. These high CPU 
times result from test problems having slow convergence near the 
optimum, which is believed to be due to many variables going to zero. 
Thus, a technique to drop variables going to zero from the LP could 
well speed up convergence. 

The ICCG algorithm does not perform very well on the test prob- 
lems. Computational results are generally inferior to those obtained 
with the MDOC algorithm. Thus, no further research into this method 
seems warranted. . 

In view of the numerical problems encountered when solving the 
least-squares problem with the normal equations approach and Cholesky 
factorization, another method that does not use square roots should be 
considered for implementation. A very promising candidate in this 
respect is Givens rotation. See Gentleman [Ref. 19] and George and 
Heath [Ref. 20]. : 
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